Instalar y cargar la librería DESeq2. Esta librería contiene una serie de funcionalidades para realizar análisis estadísticos, “aguas abajo” (down-stream) de datos ómicos analizados con SALMON e importados con tximport.
Son necesarias las librerías de “DESeq2”, “readxl”, “ggplot2”, “ggrepel” y “plotly”.
# Instalar, si no está instalado, BiocManager: cargador de dependencias de bioconductor
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
## Warning: package 'BiocManager' was built under R version 4.2.1
# Instalar, si no está instalado, DESeq2: herramienta para el análisis transcriptómico aguas abajo
if (!require("DESeq2", quietly = TRUE)) {
BiocManager::install("DESeq2")
}
## Warning: package 'DESeq2' was built under R version 4.2.2
## Warning: package 'S4Vectors' was built under R version 4.2.1
## Warning: package 'BiocGenerics' was built under R version 4.2.1
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Warning: package 'IRanges' was built under R version 4.2.1
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:grDevices':
##
## windows
## Warning: package 'GenomeInfoDb' was built under R version 4.2.2
## Warning: package 'SummarizedExperiment' was built under R version 4.2.1
## Warning: package 'MatrixGenerics' was built under R version 4.2.1
## Warning: package 'matrixStats' was built under R version 4.2.2
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Warning: package 'Biobase' was built under R version 4.2.1
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
# Instalar, si no está instalado, weitexl: paquete para guardar como excel dataframes
if (!require("writexl", quietly = TRUE)) {
install.packages("writexl")
}
## Warning: package 'writexl' was built under R version 4.2.2
# NOTA: si dice de actualizar paquetes, contestar que NO (n).
# Instalar, si no está instalado
# - readxl: sirve para leer excele
# - ggplot2: sirve para generar gráficos estáticos, muy usado en publicaciones
# - ggrepel: sirve para añadir funciones a ggplot2 que evitan que los textos dentro de un gráfico se superpongan
# - plotly: sirve para generar gráficos dinámicos, permite interactuar al usuario con la gráfica
# Normalmente, vienen instalado en muchas distribuciones de R, solo hace falta cargarlos con "library()"
Cargamos los datos clínicos de los pacientes y los datos del transcriptoma.
# --- Datos de expresion genética ---
# El objeto que contiene toda la información esta guardado como .RDS
txi <- readRDS("../import quant/txi.RDS")
# Trabajaremos con los counts
countData <- txi$counts
# No vamos a trabajar más con txi, así que como ocupa 35MB, lo elimino de la memoria
rm(txi)
# La matriz tiene, en los nombres de las muetras, un sufijo después de una barra baja, que es el "sample ID" que se utilizó para secuenciarlo en illumina. Para que no de problemas, se van a renombrar
# Elimino el sufijo "_Sxx" de los nombres de las columnas
lNewColnames <- unlist(strsplit(colnames(countData), "_"))[seq(1, 2 * ncol(countData), 2)]
# Asigno los nuevos nombres de las columnas a la matriz
colnames(countData) <- lNewColnames
# Muestro los primeros datos
head(countData)
## 103 105 112 113 131 151 165
## ENSG00000000003.15 280.690 269.872 437.320 354.990 565.889 474.571 393.377
## ENSG00000000005.6 0.000 1.000 2.000 5.000 14.000 0.000 1.000
## ENSG00000000419.14 468.499 657.489 417.199 670.997 476.536 303.273 464.651
## ENSG00000000457.14 292.789 289.795 250.316 338.842 339.367 258.275 297.864
## ENSG00000000460.17 148.557 101.415 201.240 132.638 132.400 169.704 270.774
## ENSG00000000938.13 37.000 84.008 82.001 36.000 87.000 108.000 39.180
## 17 171 177 188 201 203 210
## ENSG00000000003.15 537.758 248.298 295.417 678.890 518.671 377.485 602.191
## ENSG00000000005.6 4.000 1.000 0.000 0.000 0.000 0.000 1.000
## ENSG00000000419.14 353.991 540.465 459.604 386.604 494.334 247.871 623.557
## ENSG00000000457.14 259.594 223.485 276.667 319.040 335.244 197.801 352.266
## ENSG00000000460.17 159.832 141.252 192.590 143.592 247.313 121.909 272.588
## ENSG00000000938.13 82.001 37.000 49.000 38.000 103.999 114.000 57.000
## 218 30 37 5 63 64 70
## ENSG00000000003.15 636.489 418.424 410.440 597.728 460.811 463.123 625.235
## ENSG00000000005.6 2.000 3.000 0.000 0.000 1.000 1.000 1.000
## ENSG00000000419.14 616.128 881.789 460.925 457.024 566.649 324.236 852.521
## ENSG00000000457.14 370.288 464.448 238.741 257.118 287.889 256.463 359.007
## ENSG00000000460.17 302.957 212.938 222.348 178.380 145.989 125.019 226.644
## ENSG00000000938.13 18.000 88.000 70.999 144.001 44.000 59.000 144.590
## 77 78 80 83 88 95 98
## ENSG00000000003.15 438.495 425.433 250.667 416.714 3950.941 661.838 572.685
## ENSG00000000005.6 0.000 2.000 4.000 0.000 15.000 0.000 0.000
## ENSG00000000419.14 605.586 323.024 824.963 543.388 4038.405 825.350 526.198
## ENSG00000000457.14 313.623 232.061 315.194 242.842 2167.094 286.643 397.920
## ENSG00000000460.17 207.402 214.585 109.392 204.298 1210.723 161.989 267.191
## ENSG00000000938.13 69.000 91.030 244.000 47.001 277.017 92.343 77.000
rm(lNewColnames)
# --- Datos clínicos ---
# Cargar la librería necesaria para leer archivos excel
library("readxl")
## Warning: package 'readxl' was built under R version 4.2.1
# Leer los datos del excel
clinicalData <- read_excel("../datos clinicos/dfClinical.xlsx")
# Muestro los primeros datos
head(clinicalData)
## # A tibble: 6 × 21
## Samples Age.at.diagnosis Gender Tumor.localization Survival Grade Ki67 LOH
## <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 177 63 F Left frontal lobe 118 IV N.D. 0
## 2 80 71 F Left frontal lobe 244 IV 25 0
## 3 88 45 M Left frontal lobe 3384 II 5 N.D.
## 4 5 53 M Left frontal lobe 1239 IV N.D. 0
## 5 70 67 M Left fronto-pariet… 406 IV N.D. 0
## 6 83 52 F Left parietal lobe 22 IV 80 N.D.
## # … with 13 more variables: `IDH1.(DM)` <chr>, `IDH2.(DM)` <chr>, MGMT <chr>,
## # pTERT.C228T <chr>, pTERT.C250T <chr>, NLR <chr>, PLR <chr>, EGFRvIII <chr>,
## # fAge <chr>, fGender <chr>, fTumorHemisphere <chr>, fTumorLobule <chr>,
## # fSurvival <chr>
Para evitar errores a la hora de fusionar los datos de la cuantificación y los datos de los pacientes, se va a reordenar el dataframe de los datos clínicos de los pacientes de forma que tengan el mismo orden que las columnas de la matriz con los counts para cada gen.
# Los datos de la matriz estan ordenados como si los números fuesen carácteres, por lo que voy a transformar la variable "Sample" del dataframe de los datos clínicos a carácter, y después lo ordenaré de menor a mayor
clinicalData$Samples <- as.character(clinicalData$Samples)
clinicalData <- clinicalData[order(clinicalData$Samples),]
Elimino los datos del paciente 88, ya que por un error en el laboratorio, se secuenció con más profundidad de lectura y sus datos de supervivencia son outliers.
# Elimino la columna del paciente 88 de la matriz "countData"
countData <- countData[,colnames(countData) != "88"]
# Elimino la fila del paciente 88 en los datos clinicos
clinicalData <- clinicalData[clinicalData$Samples != "88",]
# Nos aseguramos de que hay el número correcto de muestras
ncol(countData)
## [1] 27
nrow(clinicalData)
## [1] 27
# Mostramos
head(countData)
## 103 105 112 113 131 151 165
## ENSG00000000003.15 280.690 269.872 437.320 354.990 565.889 474.571 393.377
## ENSG00000000005.6 0.000 1.000 2.000 5.000 14.000 0.000 1.000
## ENSG00000000419.14 468.499 657.489 417.199 670.997 476.536 303.273 464.651
## ENSG00000000457.14 292.789 289.795 250.316 338.842 339.367 258.275 297.864
## ENSG00000000460.17 148.557 101.415 201.240 132.638 132.400 169.704 270.774
## ENSG00000000938.13 37.000 84.008 82.001 36.000 87.000 108.000 39.180
## 17 171 177 188 201 203 210
## ENSG00000000003.15 537.758 248.298 295.417 678.890 518.671 377.485 602.191
## ENSG00000000005.6 4.000 1.000 0.000 0.000 0.000 0.000 1.000
## ENSG00000000419.14 353.991 540.465 459.604 386.604 494.334 247.871 623.557
## ENSG00000000457.14 259.594 223.485 276.667 319.040 335.244 197.801 352.266
## ENSG00000000460.17 159.832 141.252 192.590 143.592 247.313 121.909 272.588
## ENSG00000000938.13 82.001 37.000 49.000 38.000 103.999 114.000 57.000
## 218 30 37 5 63 64 70
## ENSG00000000003.15 636.489 418.424 410.440 597.728 460.811 463.123 625.235
## ENSG00000000005.6 2.000 3.000 0.000 0.000 1.000 1.000 1.000
## ENSG00000000419.14 616.128 881.789 460.925 457.024 566.649 324.236 852.521
## ENSG00000000457.14 370.288 464.448 238.741 257.118 287.889 256.463 359.007
## ENSG00000000460.17 302.957 212.938 222.348 178.380 145.989 125.019 226.644
## ENSG00000000938.13 18.000 88.000 70.999 144.001 44.000 59.000 144.590
## 77 78 80 83 95 98
## ENSG00000000003.15 438.495 425.433 250.667 416.714 661.838 572.685
## ENSG00000000005.6 0.000 2.000 4.000 0.000 0.000 0.000
## ENSG00000000419.14 605.586 323.024 824.963 543.388 825.350 526.198
## ENSG00000000457.14 313.623 232.061 315.194 242.842 286.643 397.920
## ENSG00000000460.17 207.402 214.585 109.392 204.298 161.989 267.191
## ENSG00000000938.13 69.000 91.030 244.000 47.001 92.343 77.000
head(clinicalData)
## # A tibble: 6 × 21
## Samples Age.at.diagnosis Gender Tumor.localization Survival Grade Ki67 LOH
## <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 103 63 F Right fronto-parie… 131 IV N.D. 0
## 2 105 62 M Left temporal lobe 648 IV N.D. 0
## 3 112 53 M Right occipital lo… 679 IV N.D. N.D.
## 4 113 45 M Right temporal lobe N.D. II N.D. 0
## 5 131 68 M Right temporal lobe 23 IV 60 N.D.
## 6 151 46 F Right parieto-temp… 73 IV N.D. N.D.
## # … with 13 more variables: `IDH1.(DM)` <chr>, `IDH2.(DM)` <chr>, MGMT <chr>,
## # pTERT.C228T <chr>, pTERT.C250T <chr>, NLR <chr>, PLR <chr>, EGFRvIII <chr>,
## # fAge <chr>, fGender <chr>, fTumorHemisphere <chr>, fTumorLobule <chr>,
## # fSurvival <chr>
La librería DESeq2, permite fusionar los datos de los transcritos y los datos clínicos de los pacientes. De forma que por un lado tenemos la matriz con los datos de expresion génica (counts de los transcritos) y por otro lado tenemos los datos clínicos de los pacientes que son tratados como fáctores.
@see: “https://compbiocore.github.io/deseq-workshop-1/assets/deseq_workshop_1.html”
# Cargar la librería
library("DESeq2")
# Aplicar la transformación que une counts + metadados
dds <- DESeqDataSetFromMatrix(countData = round(countData), colData = clinicalData, design = ~Grade)
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
# Mostramos la información del objeto (es un objeto especial de DESeq2)
dds
## class: DESeqDataSet
## dim: 38244 27
## metadata(1): version
## assays(1): counts
## rownames(38244): ENSG00000000003.15 ENSG00000000005.6 ...
## ENSG00000291240.1 ENSG00000291292.1
## rowData names(0):
## colnames(27): 103 105 ... 95 98
## colData names(21): Samples Age.at.diagnosis ... fTumorLobule fSurvival
Ahora vamos a aplicar un Analisis de Componentes Principales (PCA) con las funciones del paquete DESeq2.
# Normalizamos los datos
# vst = Variance Stabilizing Transformation
vsdata <- vst(dds, blind = FALSE)
# Realizamos la extracción de componentes principales
pcaData <- plotPCA(vsdata, intgroup="Grade", returnData = TRUE)
# plotPCA, por defecto, genera un dataframe con una variable llamada "group" que contiene la misma información que la variable que se ha indicado en la función. Se pude eliminar.
pcaData$group <- NULL
El objeto “pcaData” contiene los valores de las 2 primeras Componentes Principales (PCs) para los 28 pacientes.
Se ha creado una función, partiendo de la original, que calcula las 4 primeras componentes principales.
# Mostramos los datos
head(pcaData)
## PC1 PC2 Grade name
## 103 -46.78127 5.9369410 IV 103
## 105 -50.69042 6.1853577 IV 105
## 112 12.29899 -12.6832694 IV 112
## 113 -55.66857 -0.6387503 II 113
## 131 10.82648 14.7033023 IV 131
## 151 16.34288 20.5112624 IV 151
# La ubicación donde guardar el gráfico
dir.create(file.path(getwd(), "img"), showWarnings = FALSE)
# Cargamos la librería ggplot2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.2
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 4.2.1
# Extraemos el atributo del porcentaje de la varianza
# (sirve para añadir el % de la varianza en los ejes del gráfico)
percentVar <- round(100 * attr(pcaData, "percentVar"))
# Hacemos el gráfico
g <- ggplot(pcaData, aes(PC1, PC2, color=Grade, label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(paste0("PC1: ",percentVar[1],"% variance")) +
ylab(paste0("PC2: ",percentVar[2],"% variance")) +
ggtitle("Patients by glioblastoma grade") +
coord_fixed();g
# Guardamos el gráfico
ggsave("./img/PCA_01_fGrade.png", plot = g)
## Saving 7 x 5 in image
rm(percentVar, g)
La función “plotPCA” del paquete DESeq2 sola genera las Componentes Principales (CP) 1 y 2. Como Luis Miguel represento los datos de los pacientes con 3 CP, es necesario calcular las otras componentes principales.
# ----------- #
# --- PCA --- #
# ----------- #
# --- Variables --- #
# Número de genes
ntop = 500
# El fáctor
intgroup = "Grade"
# --- Extracción de PCs --- #
rv <- rowVars(assay(vsdata))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(assay(vsdata)[select, ]))
percentVar <- pca$sdev^2/sum(pca$sdev^2)
intgroup.df <- as.data.frame(colData(vsdata)[, intgroup, drop = FALSE])
# --- Dataframe con las PCs --- #
pcaData2 <- data.frame(
PC1 = pca$x[, 1],
PC2 = pca$x[, 2],
PC3 = pca$x[, 3],
PC4 = pca$x[, 4],
intgroup.df,
name = colData(dds)[,1])
# Mostramos el % de la varianza original retenida (con la varianza se estima la información)
round(sum(percentVar[1:4]) * 100, 2)
## [1] 57.31
# Eliminar variables
rm(ntop, intgroup, rv, select, pca, intgroup.df)
# Mostramos el dataframe generando con las componentes principales 1 - 4
head(pcaData2)
## PC1 PC2 PC3 PC4 Grade name
## 103 -46.78127 5.9369410 -14.2498829 8.051700 IV 103
## 105 -50.69042 6.1853577 -6.6862629 -1.714114 IV 105
## 112 12.29899 -12.6832694 -0.1276191 -6.318918 IV 112
## 113 -55.66857 -0.6387503 -6.7935896 -2.988240 II 113
## 131 10.82648 14.7033023 13.7745312 16.942380 IV 131
## 151 16.34288 20.5112624 6.2481519 -15.658440 IV 151
Por defecto, la función que realiza el PCA del paquete DESeq2 selecciona los primeros 500 genes en función de su varianza. Es decir, ordena de mayor a menor varianza las variables y selecciona las primeras 500. La poca retención de información de las componentes principales (CP) se debe en parte a que se están tomando una muestra pequeña de las variables originales (de alrededor de 37.000, solo se trabaja con 500) y porque la relación entre las variables posiblemente no sea líneal.
# ----------- #
# --- PCA --- #
# ----------- #
# --- Variables --- #
# Número de genes
ntop = 2000
# El fáctor
intgroup = "Grade"
# --- Extracción de PCs --- #
rv <- rowVars(assay(vsdata))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(assay(vsdata)[select, ]))
percentVar3 <- pca$sdev^2/sum(pca$sdev^2)
intgroup.df <- as.data.frame(colData(vsdata)[, intgroup, drop = FALSE])
# --- Dataframe con las PCs --- #
pcaData3 <- data.frame(
PC1 = pca$x[, 1],
PC2 = pca$x[, 2],
PC3 = pca$x[, 3],
PC4 = pca$x[, 4],
intgroup.df,
name = colData(dds)[,1])
# Mostramos el % de la varianza original retenida (con la varianza se estima la información)
round(sum(percentVar3[1:4]) * 100, 2)
## [1] 57
# Eliminar variables
rm(ntop, intgroup, rv, select, pca, intgroup.df)
# Mostramos el dataframe generando con las componentes principales 1 - 4
head(pcaData3)
## PC1 PC2 PC3 PC4 Grade name
## 103 -69.46740 -10.39328 2.69085208 -21.063693 IV 103
## 105 -75.02701 -16.21282 -0.07726929 -17.733730 IV 105
## 112 17.23165 11.66434 -1.66402445 6.250943 IV 112
## 113 -84.04938 -3.62071 -4.26842827 -15.594621 II 113
## 131 14.63952 -14.14128 37.05831890 12.940433 IV 131
## 151 26.65543 -37.64575 -20.23258993 20.838870 IV 151
# Guardamos los datos
# write_csv(pcaData3, "./PCA_03.csv")
# Cargamos la librería ggplot2
library(ggplot2)
library(ggrepel)
# Hacemos el gráfico
g <- ggplot(pcaData3, aes(PC1, PC2, color=Grade, label = row.names(pcaData3))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(paste0("PC1: ",round(percentVar3[1] * 100, 2),"% variance")) +
ylab(paste0("PC2: ",round(percentVar3[2] * 100, 2),"% variance")) +
ggtitle("Patients by glioblastoma grade") +
coord_fixed();g
# Guardamos el gráfico
ggsave("./img/PCA_03_fGrade.png", plot = g)
## Saving 7 x 5 in image
rm(percentVar, g)
Se utilizará un paquete distinto a ggplot2 para realizar gráficas: plotly. Plotyl es un paquete utilizado para generar gráficos dinámicos, los cuales permiten al usuario interactuar con ellos. Para la visualización de estos datos, es una mejor opción que ggplot2.
# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)
## Warning: package 'plotly' was built under R version 4.2.1
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# Generamos la figura
fig <- plot_ly(
pcaData2,
x = ~PC1,
y = ~PC2,
z = ~PC3
) %>% add_markers(color = ~Grade); fig
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
# Generamos la figura, con la opción de que cada punto tenga indicado el nombre de la muestra.
fig_named <- fig %>% add_trace(text = row.names(pcaData2), hoverinfo = 'text', marker = list(color='green')); fig_named
## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
rm(fig, fig_named)
Al dataframe que se ha generado con la función del paquete DESeq2, se le van a añadir los factores que se encuentran en el dataframe de los datos clínicos de los pacientes.
# --- PCA Data 1 --- #
# Factor edad
pcaData$fAge <- clinicalData$fAge
# Factor género
pcaData$fGender <- clinicalData$fGender
# Factor hemisferio donde aparecio el tumor
pcaData$fTumorHemisphere <- clinicalData$fTumorHemisphere
# Factor lóbulo donde aparecio el tumor
pcaData$fTumorLobule <- clinicalData$fTumorLobule
# Factor supervivencia
pcaData$fSurvival <- clinicalData$fSurvival
# --- PCA Data 2 --- #
# Factor edad
pcaData2$fAge <- clinicalData$fAge
# Factor género
pcaData2$fGender <- clinicalData$fGender
# Factor hemisferio donde aparecio el tumor
pcaData2$fTumorHemisphere <- clinicalData$fTumorHemisphere
# Factor lóbulo donde aparecio el tumor
pcaData2$fTumorLobule <- clinicalData$fTumorLobule
# Factor supervivencia
pcaData2$fSurvival <- clinicalData$fSurvival
Utilizando los factores de los datos clínicos de los pacientes, podemos comprobar si alguno de ellos muestra un patrón.
Gráficos de tres componentes principales.
Fáctor hemispherio.
# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)
# Generamos la figura
fig <- plot_ly(
pcaData2,
x = ~PC1,
y = ~PC2,
z = ~PC3
) %>% add_markers(color = ~fTumorHemisphere); fig
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
rm(fig)
Fáctor lóbulo cerebral.
# Representación de los individuos con las 3 primeras componentes principales.
# Cargamos la librería plotly
library(plotly)
# Generamos la figura
fig <- plot_ly(
pcaData2,
x = ~PC1,
y = ~PC2,
z = ~PC3
) %>% add_markers(color = ~fTumorLobule); fig
rm(fig)
Gráficos de dos componentes principales.
# Extraemos el atributo del porcentaje de la varianza
# (sirve para añadir el % de la varianza en los ejes del gráfico)
percentVar <- round(100 * attr(pcaData, "percentVar"))
# Xlab y Ylab
sXlab <- paste0("PC1: ", percentVar[1], "% variace")
sYlab <- paste0("PC2: ", percentVar[2], "% variace")
# Factor edad
g1 <- ggplot(pcaData, aes(PC1, PC2, color=fAge, label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(sXlab) +
ylab(sYlab) +
ggtitle("Patients by age quartils") +
coord_fixed()
# Factor hemisferio cerebral
g2 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorHemisphere, label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(sXlab) +
ylab(sYlab) +
ggtitle("Patients by tumor hemisphere") +
coord_fixed()
# Factor lobulo cerebral
g3 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorLobule, label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(sXlab) +
ylab(sYlab) +
ggtitle("Patients by tumor lobule") +
coord_fixed()
# Factor lobulo + hemisferio cerebral
g4 <- ggplot(pcaData, aes(PC1, PC2, color=fTumorLobule, shape = fTumorHemisphere, label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(sXlab) +
ylab(sYlab) +
ggtitle("Patients by tumor location") +
coord_fixed()
# Factor superviviencia
g5 <- ggplot(pcaData, aes(PC1, PC2, color=as.factor(fSurvival), label = row.names(pcaData))) +
geom_point(size=3) +
geom_label_repel(size = 3) +
xlab(sXlab) +
ylab(sYlab) +
ggtitle("Patients by survival") +
coord_fixed()
# Mostramos todos los gráficos
g1;g2;g3;g4;g5
# Guardamos los gráficos
ggsave("./img/PCA_02_fAge.png", plot = g1)
## Saving 7 x 5 in image
ggsave("./img/PCA_03_fHemisphere.png", plot = g2)
## Saving 7 x 5 in image
ggsave("./img/PCA_04_fLobule.png", plot = g3)
## Saving 7 x 5 in image
ggsave("./img/PCA_05_fHemisphere_fLobule.png", plot = g4)
## Saving 7 x 5 in image
ggsave("./img/PCA_06_fSurvival.png", plot = g5)
## Saving 7 x 5 in image
rm(percentVar, sXlab, sYlab, g1, g2, g3, g4, g5)
Guardar como los dataframe con las componentes principales como exceles. También, se guarda como un archivo “.RDS”.
# --- Guardar como Excel --- #
library("writexl")
# Guardamos las componentes principales 1
write_xlsx(pcaData, "./pcaData.xlsx")
# Guardamos las componentes principales 2
write_xlsx(pcaData2, "./pcaData2.xlsx")
# --- Guardar como RDS --- #
saveRDS(pcaData, "./pcaData.rds")
saveRDS(pcaData2, "./pcaData2.rds")
Borramos todas las variables
rm(pcaData, pcaData2, dds, countData, vsdata, clinicalData)